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ABSTRACT: 

In  this  report,  service  routines  are  developed  for  linear,  quadratic 
and  cubic  finite  elements  for  two  and  three  dimensional  regions.  Also,  an 
equation  solver  of  very  large  capacity  is  developed  for  systems  of  linear 
equations  where  the  matrix  of  coefficients  is  symmetrical,  positive  def- 
inite and  tightly  banded  along  the  main  diagonal. 
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INTRODUCTION 

The  stress-strain  analysis  of  solid  propellant  rocket  motors  has 
been  the  object  of  numerous  investigations  in  recent  years.    Various 
approximate  methods  of  analysis  have  been  proposed.    The  Finite  Ele- 
ment method,  without  a  doubt,   showed  the  greatest  promise  of  providing 
the  information  needed  in  the  evaluation  of  real  motors .    In  principle 
the  method  can  provide  analyses  of  any  three  dimensional  bodies.    Any 
number  of  material  and  geometrical  discontinuities  can  be  admitted. 
The  method  also  allows  for  almost  any  type  of  static,  dynamic  and 
thermal  loadings  to  be  applied.     In  practice,  however,  this  has  proved 
to  be  a  formidable  computational  task  even  for  the  largest  digital  com- 
puters in  existence. 

An  elegant  compromise  was  obtained  with  the  development  of  axi- 
symmetric  ring  elements  (1)    .     However  this  compromise  eliminated  the 
possibility  of  solving  for  real  motors  with  star  shaped  grains.    In  those 
cases  it  has  generally  been  assumed  that  a  plane  strain  analysis  using 
elements  of  triangular  shape  would  provide  enough  information  for  de- 
sign.   The  constant  strain  finite  element  triangle  has  been  used  ad 
nauseam  in  these  studies,  in  spite  of  the  fact  that  many  other  elements 
have  been  shown  to  provide  much  better  results  for  the  same  computational 
effort  (2).    When  the  perforation  changes  shape,  the  combination  of  ring 


* 
Numbers  in  parenthesis  refer  to  the  reference  list  at  the  end. 


elements  and  constant  strain  triangle  cannot  give  reliable  results. 

Recently  the  introduction  of  isoparametric  elements  in  the  literature 
has  given  a  new  impetus  to  attempts  at  solving  real  three-dimensional 
motors  (3),  (4),  (5).    It  is  to  such  a  task  that  we  address  ourselves. 
In  this  work  we  shall  not  attempt  to  account  for  real  viscoelastic  be- 
havior, what  we  shall  present  is  a  linear  elastic  analysis  method  for 
three-dimensional  objects.    We  will  show  that  isoparametric  elements 
can  be  economically  generated  in  the  computer,  even  for  the  monstrous 
(96  x  96),  32  nodal  point  bricks.    A  massive  system  of  linear  equations 
with  a  large  half -bandwidth  normally  results  from  the  analysis  of  real- 
istic problems.    For  these  systems  of  equations  we  propose  a  new 
solver  (6) . 


Chapter  I.    Isoparametric  Elements 

In  this  chapter  we  describe  isoparametric  elements  and  the  numeri- 
cal techniques  used  to  construct  them. 
1 . 1  Advantages  of  isoparametric  elements 

The  development  of  isoparametric  elements  marked  an  important 
breakthrough  in  the  finite  element  method  (4).    Until  then,  a  great 
variety  of  elements  were  developed  to  solve  specific  types  of  problems. 
The  most  successful  elements  satisfied  all  the  requirements  to  insure 
convergence  to  an  exact  result  by  mesh  refinement.    However  the  solu- 
tion of  many  problems  required  a  combination  of  various  elements  which 
could  not  satisfy  all  the  required  conditions  for  successful  convergence, 
when  assembled  to  solve  one  problem.    An  example  of  such  a  situation 
is:    an  axisymmetric  system  where  truncated  conical  ring  elements  are 
used  to  model  a  case  around  a  rocket  motor  and  ring  elements  of  tri- 
angular cross  section  are  used  to  model  the  grain  (1).    The  usual  coni- 
cal element  admits  transverse  displacements  given  by  a  cubic  poly- 
nomial where  as  the  grain  ring  elements  will  allow  only  linear  dis- 
placement in  the  same  direction.    Obviously,  adjacent  ring  and  case 
elements  cannot  maintain  continuity  of  displacements  under  strain. 

Isoparametric  element  families  on  the  other  hand  can  be  developed 
to  yield  compatible  elements  of  any  type:    line,  surface  or  volume. 
Furthermore  these  elements  admit  curved  boundaries,  resulting  in  accu- 
rate modeling  of  complex  shapes  even  with  very  few  elements.    The 


results  obtained  in  references  (3),  (4)  and  (5)  show  that  good  results 
are  obtained  with  very  coarse  mesh. 
1 .2    Stiffness  matrix  of  a  finite  element 

The  strain  energy  of  an  element  having  a  stiffness  matrix  [K]  can 
be  written  as: 

Us=^<u.>   [K]  {u.}  (1) 

where  the  vector  {u.}  represents  all  the  nodal  coordinates  of  the  element, 

The  same  energy  can  also  be  written  from  fundamental  principles 
as  follows: 


Us=I 


v 


<cr.>   U.}  dv 
i         i 


(2) 


where  {a.}  is  a  vector  of  all  the  stress  components  and  {c.}  a  vector 
of  corresponding  strain  components  in  the  volume  V.    The  constitutive 
laws  can  be  written  in  matrix  form  as  follows: 


{a.}  =  telle.} 


(3) 


where  [ft]  is  a  symmetric  matrix  of  material  properties.    The  strain  dis- 
placement relations  are  also  written  in  matrix  form  as  follows: 


{c1}  =  M< 


u(x,  y,  z) 
v(x,  y,  z) 
w(x,  y,  z) 


(4) 


where  matrix  [(p]  is  in  general  a  rectangular  matrix  of  differential  opera- 


tors, and  the  vector  <u(x,  y,  z),  v(x,  y,  z),  w(x,  y,  z)>     represents 


the  three  components  of  displacements  in  an  appropriate  global  system 
of  reference,  (x,  y,  z).    Finally  the  displacement  functions  are  as- 
sumed to  be  of  the  following  form: 

f  \ 

u(x,  y,  z) 


^  v(x,  y,  z)}  =    [n]  {u.} 
w(x,  y,  z) 


(5) 


; 


where  [ft]  is  a  rectangular  matrix  of  carefully  selected  functions  of 
(x,  y,  z).    These  functions  must  be  selected  to  insure  compatibility 
across  common  boundaries,  to  insure  that  rigid  body  motions  will  result 
in  strainless  states  and  finally  to  admit  at  least  constant  straining 
states. 

After  substitution  of  (5),   (4)  and  (3)  into  (2)  and  using  the  definition: 

En*]  =  [<p]  [n]  (6) ' 

we  obtain  the  following  result: 


uc=i  <^>(F  Cft*]T[e][n*]dv){u.}  (?) 

s       I        i     \J  v  /      1 


The  stiffness  matrix  is  then: 

[k]=  f    [n*]T[e][n*]dv  (8) 

1 .3   The  isoparametric  concept 

In  isoparametric  elements,  the  displacements  are  obtained  in  terms 
of  shape  functions  defined  in  a  system  of  curvilinear  coordinates  particu- 
lar to  one  element. 


u(x,  y,  z)  =  N  (£,n,C  )u 

v(x,  y,  z)  =  N.(C,n,c)v.  (9) 

1  1 

w(x,  y,  z)  =  N.(?  nsc)w. 

where  N.(£,n,c)  are  the  shape  functions  and  u. ,  v, ,  w.  are  nodal 
values  of  the  displacement  components  in  the  global  system  (x,  y,  z). 
The  global  system  of  reference  is  then  related  to  the  curvilinear  co- 
ordinates by  the  same  shape  functions  as  follows: 

y(5,n,c)  =N.U,n,c    )y.  (10) 

where  x.  .  y    and  z.  are  the  nodal  coordinates  of  the  element  in  the 
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global  system  (x,  y,  z).     It  is  not  our  intention  here  to  trace  the  history 
of  such  a  concept.    The  readers  are  referred  to  references  (3)  and  (4) 
for  the  details  and  justification  of  the  concept.    The  remaining  task  is 
to  construct  the  shape  functions.    We  have  taken  these  functions  from 
reference  (4).    They  will  be  found  in  section  1.6,  1.7  and  1.8. 
1.4    The  Jacobian  Matrix 

In  equation  (6) ,  the  derivatives  of  the  shape  functions  in  terms  of 
the  coordinates  of  the  global  system  (x,  y,  z)  are  needed.    These  opera- 
tions require  the  Jacobian  of  the  transformation: 
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8N. 


9t? 

ac; 


dx 

a6 

ay. 

§x 

ay_ 

dx 

ac 

dz 

ac 

'bN> 


dx 


aN. 
i 

ay 


8N. 


az 


(ii) 


where  the  square  matrix  is  the  Jacobian  [j].    Using  equation  (10)  the 
Jacobian  matrix  is  easily  evaluated  as  follows: 


[J]    = 


dN 

aN? 

H    ' 

a4 

3N 

aN? 

a?7 

a?7 

aN 

8N2 

s 


SC 


ac 


xi  yi  zi 


X2    y2    Z2 


(12) 


When  the  Jacobian  is  known  the  required  derivatives  are  easily  ob- 
tained.   The  element  of  volume  is  then  transformed  as  follows: 

dV  =  det  (J)  d£  drj  d£  (13) 

At  this  time  everything  is  ideally  suited  for  Gaussian  numerical  inte- 
gration.   This  will  be  explained  later. 
1 .  5  Shape  functions 

As  mentioned  previously  we  take  the  displacement  functions  directly 
from  reference  (4) .    We  write  these  shape  functions  for  three  dimensional 


elements  since  they  automatically  contain  the  information  for  two  and 
one  dimensional  elements.    We  also  restrict  our  attention  to  tri-linear, 
tri-quadratic  and  tri-cubic  polynomial  expressions 
1 .  6  Linear  elements 

The  shape  functions  are: 

N .  =^(1  +  4  )(1  +7?  )U  +  C)  iVl,  2 8       (14) 

where  £     =  ££    anc*  £t  is  ±  **  similarly  for  77    and  £    . 


2  (-1,1,1) 


6  (-1,1  ,-l) 


8  (1,-1,-D  5(1,1,-1) 

Figure  1.     Linear  element  (8  nodal  points) 


1 .  7  Quadratic  elements 

See  figure  2  for  the  identification  of  the  nodal  points. 
Corner  nodes:  1,  3,  5,  7 ,   13,   15,   17  and  19 

N.  =  7  (1  +  i  Ml  +  t?  Ml  +  C   )(€    +  r?    +  r     -  2)       (15) 

10  0  0  0000 


wh 


ere  4     =  ££•  and  £.  is  ±  1,  similarly  for  77    and  £    . 

nil  OO 


Mid  side  nodes:  2,  6,   18  and  14 


n.  =  7(1  -  i2)d  +n  )d  +  c) 

1    ^  00 


(16) 


Mid  side  nodes:    4,  8,  20  and  16 

N^I*1  -r?2)(i  +  C)U  +  4  ) 

1         H  0  0 

Mid  side  nodes:  9,   10,   11  and  12 


N.  =j(l  -  t2)(l  +  4   )(1  +  T7  ) 
1       q  oo 


(17) 


(18) 


Figure  2.     Quadratic  element  (20  nodal  points) 


1 .  8  Cubic  elements 

See  figure  3  for  the  identification  of  the  nodal  points . 
Corner  nodes  1,  4,  7,   10,  21,  24,  27  and  30 

n.  =  -^j  (i  +  4  )(i  +  T?  )d  +  C  )[9(£2  +  n  +  C2) 

1       Dq  o  o  o 

Mid  side  nodes:  2,  9,  29,  22,  3,  8,  28  and  23 


-  19]  (19) 


-^-(1  -  i2)(l  +  9£  )(1  +  T|  )(1  +  C   ) 

0  0  0 


i       64 


(20) 


where 


4.  =  ±  J  ,  \  =  ±  1  and  q  =  ±  1 


Mid  side  nodes:     12,  32,  25,   5,   11,  31,  26  and  6 

N.  =Jf(l  -T]2)(l  +  9t]  Ml  +Z  .)(!  +  €) 


(21) 


where 


?7.  =  ±7  ,  4.  =  ±  1  and  C.  =  ±  1 
'1  3         1  1 


Mid  side  nodes:  13,   14,   15,   16,   17,   18,   19  and  20 

N    =-~-U  =  C2)(1  +  9C   )(1  +  4.  )(1  +  t|  ) 


where 


i       64 
Ci  =  ±  3/  4j_  =  ±1  and  tj.  =  ±1. 


(22) 


30  3  I  32  21 

Figure  3.     Cubic  element  (32  nodal  points) 

In  all  of  these  formulas  if  one  coordinate    ,  for  example  C  =  C     =  1  / 

0 

the  formulas  degenerate  to  the  corresponding  correct  formulas  for  the 
two  dimensional  element  in  plane  (£ ,  77).    The  same  is  true  for  line 
elements. 
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1 .9   Numerical  integration 

When  all  the  information  contained  in  sections  1.4  to  1.8  is 
substituted  into  equation  (8)  we  finally  obtain  an  expression  of  the 
following  form: 


1      1      1 
[K]  =  [k(£,  T?,  C)l  d^drjdC  (23) 


-1  -1--1 


In  this  last  expression  the  matrix    k(£,  77,  £)  is  composed  of  elements 
in  which  the  inverse  of  the  Jacobian  matrix  is  involved,  furthermore, 
all  the  elements  must  be  multiplied  by  the  determinant  of  the  Jacobian 
matrix.    Obviously  the  process  of  reduction  to  the  form  shown  in 
equation  (23)  is  not  trivial,  however,  a  computational  algorithm  is  not 
too  difficult  to  obtain. 

Once  the  form  of  equation  (23)  is  obtained,  Gaussian  integration 
is  straightforward.    The  number  of  Gauss  points  to  be  used  in  this 
integration  is  not  easy  to  determine,  especially  if  the  element  has 
irregular  boundaries.    For  trirectangular  elements  ,  two,  four  and  five 
Gauss  ordinates  in  each  direction  of  the  local  system  of  reference 
will  give  exact  results  for  the  linear,  quadratic  and  cubic  elements 
respectively.    For  elements  having  curved  boundaries  the  same  num- 
ber of  Gauss  ordinates  can  only  give  approximate  values  to  the  stiff- 
ness matrix.    We  experimented  with  various  numbers  of  Gauss  points 
in  those  cases  where  the  results  were  approximate,  we  observed  that 
the  dominant  eigenvalues  of  the  stiffness  matrices  were  not  sensitive 


11 


to  the  number  of  Gauss  ordinates.    We  also  observed  that  the  trace 
of  an  approximate  stiffness  matrix  was  always  smaller  than  the  trace 
of  the  exact  matrix.    By  increasing  the  number  of  Gauss  points,  the 
trace  of  the  approximate  matrix  always  approached  the  exact  value 
from  the  low  side.    This  error  may  actually  be  beneficial,  since  the 
trace  of  a  stiffness  matrix  is  known  to  be  always  greater  than  the 
trace  of  the  actual  stiffness  of  a  real  structure.    For  this  reason  we 
have  adopted  to  use  2,  4  and  5  Gauss  ordinates  in  the  generation  of 
the  stiffness  matrices.    It  will  be  necessary  to  show  with  actual 
problems  whether  or  not  this  compromise  is  valid  in  general. 
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Chapter  II.    Computer  Programs  for  the  Generation  of  Stiffness  Matrices 

In  this  chapter  we  describe  several  computer  programs  used  to 
evaluate  the  stiffness  matrices  of  linear,  quadratic  and  cubic  ele- 
ments for  two  and  three  dimensional  bodies. 
2  . 1   Stiffness  matrices  for  two  dimensional  bodies 

Appendix  one  contains  the  listing  of  a  computer  program  that  was 
used  to  evaluate  the  stiffness  matrices  of  two-dimensional  elements. 
Although  the  listing  contains  only  one  quadrature  formula  with  five 
Gauss  ordinates,  an  exhaustive  series  of  tests  was  conducted  with 
formulas  using  from  two  to  ten  Gauss  ordinates.    The  five  point  formula 
was  finally  adopted  for  practical  considerations;  the  formation  times 
for  one  element  were:    0.3,   1.0,  2.5  seconds  for  linear,  quadratic 
and  cubic  elements  respectively,  and  this  was  judged  to  be  adequate 
for  this  phase  of  our  research. 

Following  the  strategy  developed  in  section  1.9,  two  sub- 
routines are  needed  to  do  all  the  work.    The  first  subroutine  named 
QUAD  5  performs  the  numerical  integration,  the  second  named  FORMK 
assembles  the  stiffness  matrix  k(£,  77,  C)  where  (£,  77,  £)  take  their 
values  at  one  of  the  Gauss  points.  Upon  exit  from  the  subroutines,  the 
user  will  obtain  the  stiffness  matrix  of  one  of  the  elements  in  this 
family  and  also  the  transformation  matrices  between  nodal  strains 
(stresses)  and  nodal  displacements.    The  calling  statement  is  given 
in  the  appendix. 
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2  .2   Stiffness  matrices  for  three  dimensional  bodies 

Appendix  two  contains  three  separate  programs  for  the  generation 
of  the  stiffnesses  of  linear,  quadratic  and  cubic  elements  for  three 
dimensional  bodies.    The  three  programs  are  similar  in  structure  to 
the  above  program,  they  were  coded  separately  for  convenience  during 
the  testing  phase  of  this  study.     The  FORTRAN  calling  statements  are 
all  contained  in  the  appendix.    The  number  of  Gauss  point  used  was 
determined  in  section  1.9  for  these  elements. 

Comments  distributed  throughout  these  programs  should  make  the 
codes  intelligible  to  an  experienced  user  of  finite  element  codes. 
2  .3   Numerical  testing  of  stiffness  matrices 

The  subroutines  mentioned  in  the  previous  two  articles  were 
tested  under  the  control  of  another  program  not  listed  here.    Eigen- 
values and  eigenvectors  were  calculated  for  all  these  elements.    Two 
dimensional  elements  always  had  three  zero  eigenvalues  and  three 
dimensional  elements  always  had  six  zero  eigenvalues.    Static  equili- 
brium checks  were    also  applied  to  all  elements  and  equilibrium  was 
always  satisfied. 
2  .4   Hybrid  elements 

A  remarkable  property  of  isoparametric  elements  is  that  once  a 
high  order  element  is  obtained  it  can  easily  be  degenerated  to  a  lower 
order  by  simple  matrix  transformations.    For  example  let  us  consider 
the  eight  nodal  point  two-dimensional  element: 
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{F}  =  [K]  {u.} 


(24) 


where  [K]  is  the  stiffness  matrix  and  {u  }  is  the  vector  of  nodal  dis- 
placements; 

{u.}  =  (VVVVVVVVVVVVV  VVV8    9 

(25) 

see  figure  4 . 

To  degenerate  this  element  to  the  corresponding  four  nodal  point 
element,  it  is  only  necessary  to  impose  a  simple  constraint  to  the  mid 
points  of  the  element.     If  the  mid-points  are  such  that  their  nodal  dis- 
placements are  equal  to  the  average  of  the  displacements  of  their 
corresponding  corner  points  the  following  transformation  matrix  is 
readily  constructed: 


where 


Cu.}=  [THuj 


[ui}  =  <u1,vlfu3,v3,u5,v5,u7,v7>- 


(26) 
(27) 


4   O 


2 


<»  8 


Figure  4,  Hybrid  elements 
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The  matrix  [T]  in  this  case  is  simply: 


[T]  = 


1.0 
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0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

0.0 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

0.0 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

0.0 

0.0 

0.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

(28) 
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The  stiffness  matrix  for  the  four  nodal  point  element  is  simply: 

[K*]  =  [T]T  [K][T]  (29) 


If  a  five  nodal  point  element  is  desired,  for  example,  suppose  we 
would  like  to  keep  nodal  point  8  in  the  above  example,  then  two  columns 
would  be  added  to  the  matrix  [T] .    Column  9  of  [  T]  would  have  a  one 
in  row  15  and  column  10  a  one  in  row  16,  all  the  other  terms  of  rows 
15  and  16  would  be  zeros.    The  congruent  transformation  (29)  would 
give  the  stiffness  of  the  hybrid  element. 

With  such  elements  it  is  possible  to  combine  different  types  of 
elements  in  the  solution  of  one  problem  and  still  preserve  all  the 
conditions  guaranteeing  convergence  to  an  exact  result  by  mesh  size 
reduction. 
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Chapter  III.    An  Equation  Solver  of  Very  Large  Capacity 
In  this  chapter  a  very  simple  equation  solver  of  very  large  capa- 
city is  developed.    The  method  used  does  not  require  a  very  large  in- 
core  memory;  however,  a  fast  random  access  device  like  a  magnetic 
disk  drive  is  required. 
3 . 1   Introduction 

The  solution  of  very  large  systems  of  linear  equations, 

{R}  =  [K]  •   [r]  (30) 

where  [K]  is  symmetric  and  positive  definite,  is  an  activity  that  must 
be  entered  by  the  stress  analyst  when  using  the  finite  element  method. 
When  the  first  automatic  stress  analysis  programs  were  developed,  the 
Gauss-Seidel  iteration  method,  with  various  over  relaxation  schemes,, 
was  almost  universally  employed  to  solve  the  equations.    It  was  soon 
realized  that,  because  of  the  properties  of  [K],  a  direct  Gaussian 
elimination  algorithm  would  be  much  more  economical  of  computer  time. 
However,  matrix  [K] ,  must  be  in  core  during  elimination,  and  this  puts 
a  very  stringent  limitation  on  the  total  number  of  equations  that  can 
be  solved. 

For  the  purpose  of  comparing  various  methods,  we  will  presume 
that  we  have  at  our  disposition  a  machine  with  32,000  words  of  core 
storage.    Of  the  total  memory/  let  us  assume  that  approximately 
25,000  words  are  available  for  the  storage  of  [K].    For  such  a  machine, 
then,  approximately  160  equations,  at  most,  can  be  solved  by  this 
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method. 

Because  the  matrix  [K]  is  usually  sparse,  a  special  ordering  of 
the  unknowns  can  produce  a  matrix  [K]  that  is  tightly  banded  about 
the  main  diagonal.    Taking  advantage  of  this  character  of  [K],  several 
programs  were  devised  that  required  only  (N.M)  words  for  the  storage 
of  [K].     (See  Figure  5  for  the  meaning  of  N  and  M) .    The  short  program 
listed  in  Appendix  3  is  of  that  character.    If  M  is  very  small,  then 
the  number  of  equations  that  can  be  solved  can  become  fairly  large. 
For  example,  using  the  same  machine  as  above,  if  M  =  50,  then 
N  =  500. 

With  the  development  of  two  and  three  dimensional  elements, 
however,  this  improvement  proved  insufficient.    These  problems  fre- 
quently lead  to  systems  of  several  thousands  of  equations  with  a 
half  band  width  in  the  hundreds .    Equation  solvers  were  then  developed 
that  did  not  require  to  have  all  the  equations  in  core.    Professor  E. 
Wilson  of  the  University  of  California  at  Berkeley  (U .  S.)  several 
years  ago  developed  a  band  solver  that  could  be  used  with  tape  stor- 
age to  solve  large  systems.    The  method  used  required  (2M.M)  words 
of  in-core  storage  for  [K]  and  the  total  number  of  equations  that  could 
be  solved  was  only  limited  by  the  capacity  of  a  tape.    For  the  machine 
used  above,  the  maximum  value  of  M  is  approximately  110.    Subsequently, 
as  mentioned  in  Reference  (8),  several  workers  developed  algorithms 
requiring  only  (M  +  2)(M  +  3)/2  words  of  core  storage.    For  the  machine 


19 


used  above,  the  maximum  value  of  M  is  approximately  220. 

With  large  systems  of  equations,  another  difficulty  becomes 
evident.     Round  off  starts  to  play  a  very  important  role.    An  obvious 
remedy  is  to  convert  everything  to  double  precision.    However,  this 
roughly  doubles  the  core  storage  requirement.    Round  off  can  some- 
times be  eliminated  by  using  the  same  programs  to  calculate  corrections 
to  the  answer  by  an  iterative  technique.    However,  such  corrections 
are  not  always  valid  when  the  matrix  [K]  is  locally  ill-conditioned. 
For    this  reason,  many  workers  prefer  to  use  double  precision  all  the 
way.    The  block  solver  presented  here  uses  double  precision  every- 
where. 

As  mentioned  in  References  (9)  and  (10),  no  method  has  been 
found  that  is  faster  and  more  accurate  than  Gaussian  elimination  for 
a  system  in  which  the  matrix  of  coefficients  is  dense.    For  dense 
system  banded  along  the  main  diagonal,  it  is  reasonable  to  conjecture 
that  the  same  conclusion  will  hold  if  the  algorithm  is  modified  to 
eliminate  all  operations  involving  the  off  diagonal  zeros.     Improve- 
ments can  only  be  accomplished  by  taking  advantage  of  the  sparsity 
within  the  band,  as  shown  in  References  (11)  and  (12).    These  improve- 
ments, however,  are  obtained  at  the  cost  of  a  more  elaborate  code  to 
recognize  the  particular  type  of  sparsity  that  can  lead  to  important 
reductions  of  computations.    The  economies  of  storage,  on  the  other 
hand,  are  not  sufficient  to  allow  for  the  solution  of  very  large  systems. 
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Frontal  solvers,  as  described  in  References  (13)  and  (8),  also  produce 
savings  of  storage  and  execution  time,  but  again  prohibit  the  solution 
of  very  large  systems. 

Our  primary  concern  here  was  to  design  a  solver  that  would  permit 
the  solution  of  very  large  systems  even  with  a  small  computer.    The 
block  solver  proposed  in  this  paper  does  not  pose  any  limits  on  the 
half  band  width  nor  on  the  total  number  of  equations,  as  long  as 

storage  is  available  on  a  fast  random  access  external  device.    The 

* 
program  contained  in  Appendix  4  requires  5800  bytes     of  storage  for 

the  object  codes  and  a  working  space  depending  on  the  block  size 

selected;  for  block  sizes  of  (25,  50,   100)  the  storage  requirements 

are  (15600,  61200,  242400)  bytes  respectively.    As  mentioned  above, 

a  fast  random  access  facility  is  required.    In  our  computer  center, 

we  used  IMB/2314  disk  packs.    On  such  a  device,  55000  records, 

400  bytes  long,  can  be  stored.    This  is  enough  for  2,750,000  double 

words.    For  example,  5000  equations  with  a  half  band  width  of  451 

would  fit  with  a  complete  load  vector  on  one  disk. 


M  is  the  half  band  width 

N  is  the  total  number,  of 
equations 


Figure  5.  Banded  Character  of  K 


* 


A  byte  is  a  word  segment  of  8  bits  of  information. 
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3  .2   Elimination  by  blocks 

Consider  the  system  of  equations  shown  below: 


'r.^ 


R, 


R, 


RM     = 


R, 


v\ 


Kll       K12 


K12      K22 

T  T 

K13       K23 


T  T 

K1M     K2M 
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13 


K23      K24 


K33      K34 
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l-M+1,1 
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1-1,1 
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1M     ' 


^2,M+1     * 


K 


35 
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M+2 


^-^M^M   \l,M+l   '  H/L, 


2M 
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II 


\ 


1+1 


N-M+1,N 


\ 


I+M-] 


V 


rM> 


W 


N      NN 
(31) 


r 


The  first  equation  is  solved  for  r    to  get: 

rl  =  K11  (R1  "K12  r2  "K13  r3 


■kimV     (32) 


After  substitution  in  equations  2  to  M,  we  obtain  the  following  results: 


-  T       -1 

R2  "  R2  "K12  Kll  Rl 


(33) 


-  T         -1 

RM"RM  'K1MK11  Rl 


(34) 


for  the  load  vectors.    All  the  coefficients  in  rows  2  to  M  from  Column  2 
to  M  are  reduced  by  similar  operations;  a  typical  example  for  term  K    .  is: 
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—  T       -1 

K34"K34  "K13  Kll  K14  (35) 

Then,  Equation  2  is  solved  for  r    and  substituted  in  Equation  3  to  M  +  1 

and  the  entire  process  is  repeated  until  we  get  only  the  last  equation: 


RN  "  KNN  FN  (36) 

This  last  equation  can  now  be  solved  for  r     and  a  simple  process  of 
back  substitution  in  reverse  order  gives  all  the  unknowns. 

The  elimination  process  described  above  preserves  the  symmetry 
of  [K].    This  process  is  valid,  whether  the  elements  of  [K]  are 
individual  coefficients  or  square  submatrices  of  coefficients.    From  now 
on,  we  presume  that  the  [K    ]  are  submatrices  that  we  will  call  blocks. 
A  careful  examination  of  the  elimination  steps  reveals  that  no  more 
than  three  different  blocks  of  coefficients  are  needed.    With  this  in 
mind,  the  rest  of  the  task  is  trivial.    The  simple  solver  contained  in 
Appendix  3  was  modified  to  do  the  elimination  by  blocks.    Appendix  4 
contains  the  listings  of  the  four  subroutines  needed. 
3.3   Listings 

Subroutine  SOLVE  is  self-explanatory  and  needs  no  elaboration, 
except  as  to  the  manner  in  which  the  equations  must  be  stored.    The 
blocks  are  stored  sequentially  by  rows,  as  shown  below.    At  the  end 
of  each  row,  a  block  of  load  vectors  is  present: 
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KK10KK        K.  Rn 

11      12      13     14  lm  1 


K22  K23  K24  K25 K2  m  R2 


Kn  ki,i+i     ■■;..-■ RI  <37) 


K.        .  K    .         R      _ 

n-l,n-l    n-l,n  n-1 


K       R 

nn  n 


The  symbols  m  and  n  are  the  number  of  blocks  per  row  and  column, 
respectively. 

This  method  of  storage  leaves  a  few  blocks  of  unused  storage. 
However,  the  last  row  of  blocks  is  needed  for  temporary  storage  during 
the  elimination.    Three  service  subroutines  are  needed  in  SOLVE.    The 
multiplication  and  inversion  operations  are  standard  and  need  no  words 
of  explanation.    The  subroutine  WRDISK/RDDISK,  however,  is  hardware 
dependent  and  uses  FORTRAN  commands  particular  to  IBM/360  installa- 
tions (Reference  14).    This  feature,  although  undesirable,  is  not  a  very 
serious  drawback  because  every  installation  using  fast  random  access 
devices  has  specialized  commands  to  use  these  facilities.    For  other 
installations,  this  subroutine  would  have  to  be  rewritten  to  conform 
with  the  available  hardware.    The  listings  are  self-explanatory. 
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3.4  Timing 

The  use  of  the  program  listed  in  Appendix  4  poses  a  serious  practi- 
cal problem.    The  solution  of  large  systems  of  equations  requires  a 
large  amount  of  computer  time.    The  most  important  factors  affecting 
the  required  time  in  the  system  proposed  are:    the  number  of  calls  to 

the  WRDISK/RDDISK  routine,  and  the  size  of  the  blocks  (N  ) .    In  the 

s 

program,  the  number  of  calls  to  WRDISK/RDDISK  has  been  reduced  as 
much  as  we  could.    In  actual  tests,  the  object  codes  were  optimized, 
with  respect  to  execution  time,  by  the  compilers  that  produced  these 
object  codes.    We  used  the  H  compiler  of  release  18  of  the  IBM/360 
operating  system.    A  system  of  500  equations  with  a  half  band  width 
of  151  was  solved  with  three  different  block  sizes,   (25,  50,   100). 
The  execution  times  were  (510,  402,  297)  seconds,  respectively.    On 
the  basis  of  this  test,  it  seems  that  the  block-size  should  be  as  large 
as  the  computer  will  allow,  in  order  to  reduce  the  execution  time.    For 
a  block-size  of  50,  the  following  empirical  formula  gives  an  estimate 

of  the  required  time: 

2  2  2  2 

T=  1.857  nm     +  0. 102  n  m  -  0.322  n    +  4. 682  nm  -  5.230  m      (38) 

The  formula  was  obtained  with  the  tests  shown  in  Table  1.    For  the 

first  test  n  =  4,  m  =  3,  it  was  possible  to  obtain  an  in-core  solution 

with  the  subroutine  shown  in  Appendix  3.    The  execution  time  was  61.3 

seconds  for  the  in-core  solution,  as  compared  to  75.7  seconds  for  the 

block  solver.    This  is  indeed  a  very  small  penalty  if  we  also  compare  the 
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storage  requirements  of  the  two  solvers;  196K  bytes  for  the  in-core 
solver  and  122K  bytes  for  the  block  solver. 

In  a  recent  study,  Reference  (15),  various  equation  solvers  were 
compared.    Using  a  sliding  block  technique  similar  to  that  mentioned 
in  Reference  (8),  a  system  of  1080  equations  with  a  half-band  width 
of  357  was  solved  in  1365.1  seconds.    The  machine  used  was  a  CDC 
6600  and  the  method  used  required  73728  words  of  core  storage. 
Test  Number  5  in  Table  1  shows  that  1000  equations  with  a  half-band 
width  of  351  required  2989.8  seconds  and  used  152  50  double  words  of 
core  storage.    Again,  this  is  a  comparison  that  seems  very  favorable 
to  the  block  solver,  if  we  keep  in  mind  that  the  CDC6600  is  a  much 
faster  machine  than  the  IBM/360/67.    Average  execution  times  for 
division,  multiplication  and  addition  are:  (2.9,   1.0,  0.4)  micro- 
seconds for  the  CDC  machine  and  (13.9,  10.6,  8.0)  for  the  IBM  machine. 
3 . 5  Accuracy 

The  most  critical  operation,  with  respect  to  round-off,  is  the 
inversion  of  the  diagonal  blocks.    If  the  block  is  singular,  the  program 
will  print  a  diagnostic  to  this  effect.     However,  this  is  not  likely  to 
happen  very  often.     A  much  more  serious  difficulty  comes  from  blocks 
that  are  nearly  singular  and  this  must  also  be  detected.    For  this 

reason,  the  inversion  routine  will  calculate  an  estimate  of  the  round-off 

-12 
value  of  zero  as  the  trace  of  the  block  multiplied  by  10      "  and  report 

that  the  block  is  nearly  singular,  if  any  diagonal  term  becomes  smaller 
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TABLE  1       TIMING  TESTS 


n  m  N 

(Number  of  (Number  of  (Number  of 

blocks  per  blocks  per  equations) 
column)                row) 


M  T 

(Half-band    (Time  in  seconds) 
width) 


4 

3 

200 

101 

75.7 

100 

3 

5000 

101 

2977.4 

10 

5 

500 

201 

586.5 

50 

5 

2500 

201 

3834.8 

20 

8 

1000 

351 

2989.8 

The  blocksize  used  was  50  and  the  total  storage  requirement  of  the 
test  program  was  122K  bytes  or  152  50  double  words. 
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than  this  value  during  inversion.     Besides  these  diagnostics,  the 
program  will  also  print  a  condition  number  for  each  diagonal  block. 
This  condition  number  is  obtained  as  the  product  of  the  maximum  column 
sum  of  the  block  before  and  after  inversion,  following  Wilkinson's 
perturbation  theory.    With  this  condition  number,  an  estimate  of  the 
relative  error  in  the  answer  can  be  obtained  from: 

E(%)   £  c  •  ||| '  •   100  (39) 

where  c  is  the  maximum  value  of  all  condition  numbers  of  the  system, 
||p||  is  the  sum     of  the  absolute  values  of  the  estimated  errors  in  the 
loads  and  ||r||  is  the  sum  of  the  absolute  values  of  the  loads. 

3 . 6  Modifications 

The  program  is  presented  in  its  simplest  form;  no  facilities  for    • 
iteration  and  multiple-load  vectors  are  provided  for.    Users  with  particu- 
lar needs  should  have  no  difficulty  modifying  the  program  listed,  in 
order  to  take  advantage  of  hardware  available  in  their  computing 
facilities. 

3. 7  Conclusion 

The  program  presented  here  provides  a  very  simple  answer  to  a 
most  difficult  practical  difficulty.    When  a  problem  is  too  big  for  a 
band  or  a  frontal  solver,  the  block  solver  will  still  give  a  solution, 
even  with  a  very  small  computer,  provided  one  is  ready  to  pay  a 
penalty  in  execution  time . 
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Chapter  IV.    Conclusions 

The  publication  of  this  report  was  delayed  in  order  to  verify  the 
integrity  of  the  service  routines  that  it  contains.    This  has  now  been 
completed  and  it  can  be  reported  that  the  routines  perform  as  expected. 

4.1  Stress-strain  analysis  codes 

Several  computer  codes  for  mesh  generation  and  consistent  load 
vectors  computation /  as  well  as  two  basic  stress  strain  programs  for 
two  and  three  dimensional  problems,  have  been  completed  and  tested 
(16),   (17)  and  (18).    From  that  work  it  may  be  concluded  that  these  ele- 
ments will  shed  significant  new  light  on  the  stress  analysis  of  rocket 
motors . 

4.2  Reduced  Gaussian  Integration 

In  both  reference  (17)  and  (18),  it  is  reported  that  the  Gaussian 
integration  mentioned  in  article  1.9  leads  to  improved  results  when 
only  2  and  3  Gauss  ordinates  are  uses,  respectively,  in  the  generation 
of  the  20  and  32  nodal  point  elements.    This  reduced  integration  also 
leads  to  significant  economies  in  the  total  solution  times. 

4. 3  Surface  stresses 

In  some  problems  it  was  observed  that  the  stress  and  strain  values 
obtained  at  the  nodes  were  in  error  on  the  external  surfaces  of  the  ob- 
jects analysed.    Internal  nodes,  however,  gave  the  best   results.    This 
matter  deserves  further  investigations. 

4.4  Future  work 

The  program  developed  during  the  course  of  this  investigation  will 
be  refined  and  released  as   soon  as  they  are  ready.    In  the  meantime, 
potential  users  can  refer  to  reference  (17)  for  a  first  release  of  the  three- 
dimensional  programs. 
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Appendix  I.    Code  for  the  generation  of  two-dimensional  elements 
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Appendix  2  .    Codes  for  the  generation  of  three-dimensional  elements 
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Appendix  3 .    Listing  for  an  in-core  band  solver 
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Appendix  4 .    Listing  for  the  Block  Solver 
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